Generic Tracking of Multiple Apparent Horizons with Level Flow 
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We report the development of the first apparent horizon locator capable of finding multiple apparent 
horizons in a "generic" numerical black hole spacetime. We use a level-flow method which, starting 
from a single arbitrary initial trial surface, can undergo topology changes as it flows towards disjoint 
apparent horizons if they are present. The level flow method has two advantages: 1) The solution 
is independent of changes in the initial guess and 2) The solution can have multiple components. 
We illustrate our method of locating apparent horizons by tracking horizon components in a short 
Kerr-Schild binary black hole grazing collision. 



I. INTRODUCTION 

Our goal is to investigate the strong field regime of gen- 
eral relativity. In particular we wish to focus on the study 
of coalescing black hole binaries. Over the last three 
decades since the pioneering work of Cadez, Smarr, Ep- 
pley and others, advances in computing technology, nu- 
merical algorithms and techniques and our understand- 
ing of the underlying physics have advanced to a point 
where we are able to carry out simulations of binary black 
hole collisions in 3+1 dimensions. One of the outcomes 
of such simulations will be an understanding of the un- 
derlying physics of the problem; and, therefore, a pre- 
diction and understanding of the gravitational radiation 
content. A detailed knowledge of how the resultant grav- 
itational waveforms relate to the physical parameters of 
the binaries that produce them will be of importance to 
gravitational wave observatories (such as LIGO, VIRGO, 
TAMA, GEO600) now under construction around the 
world. With the building of these observatories, we stand 
at the epoch of the first direct observations of astrophys- 
ical sources that involve strong field general relativity. 

The orbit and merger of two black holes is one candi- 
date source for ground based detection of gravitational 
waveforms. This is of great interest to the relativity com- 
munity. The binary black hole problem is a two-body 
problem in general relativity. It is a stringent dynam- 
ical test of the theory. However, studying spacetimes 
containing multiple black holes involves solving the Ein- 
stein equations, a complex system of non-linear, dynamic, 
elliptic-hyperbolic equations intractable in closed form. 
The intractability of the problem has led to the develop- 
ment of numerical codes capable of solving the Einstein 
equations. 

Our approach to numerically solving the Einstein equa- 
tions involves reformulating them as an initial value prob- 
lem. In this 3+1 formulation Q, spacelike hypersurfaces 
parameterized in time foliate the spacetime. The result- 
ing equations are coupled elliptic and hyperbolic differ- 
ential equations of the 3-metric, fc, and the extrinsic 



curvature, Kij. The initial value problem is solved by 
specifying a hypersurface at an instant of time, say t = 0, 
and evolving to the next hypersurface at time t = St with 
the evolution equations to obtain and at the next 
time t — St. 

One vital issue in numerically solving the Einstein 
equation describing spacetimes containing black holes is 
handling the physical singularities. As one approaches 
the singularity, the values of the fields being computed 
approach infinity; therefore, a region containing the sin- 
gularity must be avoided to keep the computation from 
halting. Excision techniques are promising in avoiding 
the singularity. Excising the singularity involves locat- 
ing a region interior to the event horizon containing the 
singularity on each evolving hypersurface. This region is 
then "masked" from the computation. The derivatives at 
the boundary between the masked region and the compu- 
tational domain are handled using causal differencing, a 
differencing scheme |?) that respects the causal structure 
of the spacetime. 

In deciding where in space to excise we use the appar- 
ent horizon as opposed to the event horizon. By its very 
nature, the event horizon is a global construct depending 
on the entire spacetime. The apparent horizon, a local, 
i.e. spacelike 2-surface is more suitable. Following a sug- 
gestion by Unruh ||, the apparent horizon is used to 
define the excised region to be masked at each time dur- 
ing the evolution. Apparent horizons are defined locally 
for each time, and exist at, or inside of, the event horizon. 
In some spacetimes, choices of foliation may lead to the 
absence of an apparent horizon. When the discussions 
in this refer to a hypersurface, it is assumed an apparent 
horizon exists on that hypersurface. 

Recent three-dimensional horizon locator codes are ca- 
pable of finding the location of an apparent horizon 
in generic single black hole spacetimes and two 

pO| , pT[ are capable of finding disjoint multiple apparent 
horizons in the special case of conformally flat binary, 
time-symmetric black hole spacetimes. Multiple appar- 
ent horizon finding algorithms will be necessary in sim- 
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ulations of generic binary black hole spacetimes. The 
method presented in this paper, called the level flow 
method, is capable of detecting multiple apparent hori- 
zons in generic spacetimes. The level flow algorithm has 
two advantages: 1) It is independent of a good initial 
guess and 2) It is capable of following the surface through 
a change in topology. In level flow, the apparent horizon 
equation is reformulated as a parabolic equation and a 
set of surfaces are flowed with speeds dependent on the 
expansion of the outgoing null vectors normal to each 
surface. 

The purpose behind developing the level flow method 
of tracking apparent horizons is to have a method capa- 
ble of detecting multiple apparent horizons on any given 
hypersurface without a good guess. Specifically, we want 
a tracker that can detect the transition from a double to 
a single apparent horizon in single time step without a 
prior knowledge of the transition. 

In the rest of this paper we discuss apparent horizons 
i n ge neral and current 3D work in § [I7| and § In 



IV we describe the level flow method in detail and give 



a brief description of the numerical method involved in 
solving the apparent horizon is in § Q|. We demonstrate 
the use of the level flow tracker on mode l data and a 
binary black hole grazing collision in § and § VI][ 



II. APPARENT HORIZONS 




The expansion of the outgoing normals, V a fc a , can be 
rewritten in terms of quantities defined on the hypersur- 
face: 



D a i 



K + s a s b K ab 



(3) 



D a denotes covariant differentiation with respect to g a b, 
is the extrinsic curvature, and K is g ab K a },. In fact, 
there is a level set of surfaces in E parameterized by k. 
Each surface in the level set is defined by the constant 
expansion of its null vectors such that 

K = c„, (4) 

where c„ are constants labeled by the positive integer n. 
Marginally trapped surfaces arc members of this set for 



k = 0. 



(5) 



Eqn. (|^) is called the apparent horizon equation since the 
apparent horizon is the outermost surface with k = in 
E. 

The S 2 topology of the apparent horizon naturally 
lends itself to characterization via spherical coordinates. 
The function, 



ip = r-h(6,<f>) 



(6) 



is a level set of 2-spheres in E, where h(9, <f>) is called the 
apparent horizon shape function. A marginally trapped 
surface has %b = 0. The spacelike normals to S are defined 
from eqn. (^|) such that 



= g ij d j ^/y/g kl d k i>di^ 



(7) 



is the spacelike vector field at every point of S. 
The apparent horizon equation in spherical coordinates 
(h(9, 4>), 9, <j)) is a 2-dimensional problem in 9 and <j>. 



FIG. 1. Representation of a 2-sphere embedded in a hyper- 
surface, E. 

M is the spacetime with metric A g a b foliated by hyper- 
surfaces E parameterized by t with 3-metric g ab . Let S 
be a surface with S 2 topology on E. The apparent hori- 
zon is the outermost marginally trapped surface in E, a 
surface with zero expansion. The zero expansion of the 
surface is defined in terms of outgoing null vectors to 5, 
k a , such that k a have zero divergence 

V a k a = 0, (I) 

where V a is the covariant derivative associated with 4 g a b- 
Referring to fig. ([!]), k a is denned in terms of the space- 
like normals to S, s a , and the future directed timelike 
normals, n a , such that: 

k a = s a +n a . (2) 



III. CURRENT 3-D METHODS 

The approaches to solving the apparent horizon equa- 
tion on a three-dimensional hypersurface can be ad- 
dressed roughly in two categories: methods that solve 
the apparent horizon equation directly and methods that 
solve it by first recasting it as a parabolic equation. This 
paper does not address spherical and axi-symmetric ap- 
proaches. One of the first three-dimensional apparent 
horizon solvers was published by Nakamura, Kojima, and 
Oohara ||. They directly solve the apparent horizon 
equation by using spherical harmonic basis functions to 
expand the apparent horizon shape function, h(9, (f>) into 

h(6,cf>) = £ Yl a lrn Y lrn {9,4,). (8) 
;=0 m=— I 

This method is called the pseudospectral method. A fi- 
nite number of the coefficients, {a; m } parameterize the 
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horizon shape function, and the maximum l max depends 
on the computation and deviations from a sphere. The 
apparent horizon equation can then be solved by writing 
it as 



I K(ai m ) || = 0, 



(9) 



and using functional integration routines to find the co- 
efficients ai m . Others have used similar methods 

In another approach to direct solutions of the appar- 
ent horizon equation is to treat it as a boundary value 
problem. One notes that a discretization of this equa- 
tion leads to a system of algebraic equations which can 
then be solved via Newton's method. Thornburg [Q dis- 
cusses applications of Newton's method to this problem 
in general and shows results in axisymmetry. Huq |l3| ] 
has implemented a similar algorithm based on Newton's 
method that utilizes Cartesian coordinates to difference 
the equations. 

Tod |l4j first suggested the use of curvature flow in 
solving the apparent horizon equation by recasting it as 
a parabolic equation. Bernstein Jl5[ implemented Tod's 
suggestion in axisymmetry. Gundlach introduced a 
fast flow method which combines the ideas of the flow 
method with the pseudospectral method. Pasch |ll[] and 
Diener JIo[| implemented a similar method, a level-set 
method, in three-dimensions and found discrete apparent 
horizons in multi-black- hole spacetimes; however these 
spacetimes were confined to be conformally flat and time- 
symmetric. 

Each of the approaches briefly described above, solving 
the apparent horizon equation directly or solving it via 
a parabolic equation, has its advantages. Direct meth- 
ods tend to be faster while flow methods do not rely on 
"good" initial guesses. However, none of these methods 
are applicable to the generic, multi-black-hole problem. 
Herein lies the motivation behind the level flow method. 
The level flow method is the only method designed to 
locate discrete apparent horizons in generic spacetimes 
containing one or two discrete horizons. 



IV. LEVEL FLOW METHOD 

A. Curvature Flow 

The level flow method is a hybrid flow/level— set 
method. The previous section mentioned the flow 
method, this section gives more detail on the flow method 
which is the foundation of the level flow method. The 
flow approach, as suggested by Tod, is to rewrite the 
apparent horizon equation as the speed function in a 
parabolic equation. In the case of a time symmetric 
hypersurface, in which K a b = 0, the apparent horizon 
equation reduces to the condition for a minimal surface, 
D a s a = 0. In this case, the surface S is at a local ex- 
ternum of the area. The starting surface, S(X = 0), is 



parameterized by coordinates x a and evolved in terms of 
a parameter A. The equation of motion is 



dx a 



= -Hs a 



(10) 



where dx a /d\ is a vector field, and H is the mean curva- 
ture, which is the trace of extrinsic curvature associated 
with embedding S in S given by 



H = D a s a . 



(11) 



Eqn. (jl(]) is the gradient flow for the area functional. The 
area decreases monotonically with increasing A. Grayson 
Jl6| has shown that a surface deforming under its gra- 
dient field (Eqn.©) will evolve to a stable minimum 
surface (surface is local minimum of area) if there is one, 
otherwise to a point. 

In numerical relativity, we are interested in the generic 
case, with K a \, ^ 0, for which the marginally trapped 
surfaces differ from minimal surfaces, the surfaces are not 
extrema of the area. However, Tod suggests an equation 
similar to Eqn.(|To|) as a curvature flow: 



dx a 



F(K)s a 



(12) 



using F(k) = — k where k — £>aS° + s a s b K a i, — K as 
in eqn.(||). We have found eqn. ( |l2| ) to be a successful 
practical implementation of the flow method. 



B. Level Flow 

Eqn. ( ^2[ ) gives us an initial value problem. Given in- 
formation about the system at some initial A, eqn. ( p^ ) 
will describe the system for all its future propagation in 
A. Directly solving eqn. ( |l2"| ) will lead to a successful de- 
tection of single apparent horizons; however, solving it 
directly does not ensure correct handling of a topology 
change which is necessary for detection of multiple appar- 
ent horizons. By combining the flow method with a level- 
set idea however, this topology change can be effected 
and multiple apparent horizons can be tracked starting 
from a single guess surface. 

First eqn. ( |l2|) is recast from an equation governing 
the motion of the coordinates parameterizing S, namely 
x a , to an equation governing the motion of the surface 
tjj. Noting that ip is parameterized by A, 



dip dx a dip 



(13) 



d\ d\ dx a 
by the chain rule, and multiplying eqn.(|l2|) by gives 
dip 



OX 



F(n)s a dip/dx a . 



(14) 



Using 



3 



and 



IIW|| 



dib dtb , 
— — 9 ab , 



dx a dx b • 



the test surface's flow is given by: 



dip 



= F(k)\\VtP\\. 



(15) 



(16) 



(17) 



is a reformulation of eqn. (|lj) and will flow the 
surface, ip, to a marginally trapped surface at ip = when 

K = 0. 

The strength of the flow methods is their ability to lo- 
cate a surface with k — given any non-pathological 
initial surface. For example, the apparent horizon in 
a spherically symmetric spacetime (Schwarzschild) was 
found by flowing an initial surface shaped as a leaf, see 
Fig. (||). This ability is especially important when track- 
ing horizons during evolutions of binary black hole space- 
times. In this case, finding apparent horizons for two dis- 
crete apparent horizons in each involves flowing two 
initial guesses, one for each horizon jl7|]. On £(t = 0), the 
location of the apparent horizons may be known; how- 
ever, as the black holes accelerate the task of guessing 
the locations of the two horizons becomes more difficult. 
Further, the two horizons merge into a single horizon at a 
single instant of time, rendering a good initial guess diffi- 
cult. Some way of determining when two horizons merge 
into a single horizon is necessary. The level flow method 
takes care of these issues by not requiring a good ini- 
tial guess (V>(A = 0)) and by detecting multiple apparent 
horizons from a single guess (^(A = 0)). 





FIG. 2. Example of a surface undergoing flow on a spher- 
ically symmetric hypersurface. The figure labeled (a) is 
ip(X — 0), the initial guess. Figure (b) is ip(X) representing 
the final solution, the apparent horizon. 

The level flow method differs from the flow method in 
the specification of the speed function, F(k). F(k) = 
determines when the propagation of the surface stops. 
The flow method is in the form of eqn. (|l2|), in which 
F(k) = — k. A good choice since F(k) = —n = indi- 
cates a marginally trapped surface; but this choice will 
not flow 4> though a fission. In general the scheme fails 
as the surface pinches due to ill-defined normals at the 
surface. The level flow method alleviates this problem by 
looking for indications that the surface topology is about 
to change before the pinching occurs. (Another method 
which was introduced by Sethian and Osher |l8[ ] for non- 
relativistic problems is to flow a higher dimensional sur- 
face in which ip is embedded. This higher dimensional 
surface does not fission. This has only been implemented 
in a time-symmetric spacetime |flljl and requires more 
computational power due to the extra dimension in the 
problem.) 

The level flow method flows a set of two-dimensional 
surfaces in £ parameterized by k. We call the set of 
surfaces a level set and label the set S(c n ). Each surface 
has a constant value of k = c n everywhere on it. The set 
of surfaces is defined by varying c„ as the flow progresses 



± Ac, 



(18) 



where (+) indicates outward flow, (— ) inward flow, and 
Ac oc || k || 2. Each surface obeys the equation of motion 
given in eqn. ( \n\) with F(tt) defined to flow to multiple 
surfaces. We choose two options for the speed function: 



F(k) = k — c n 

F(k) — (k— c n ) arctan 2 (- — — ). 



(19) 
(20) 



As k — c n — * 0, both functions are solving for a partic- 
ular surface in the level set, S{c n ). The second func- 
tion, eqn.(|2fj|), behaves similarly to the first but allows 
for larger time steps near a fissioning surface because it 
moves points further from the k — c n = surface faster 
than the points closer to this surface. 

Eqn. (|17J) is an initial value problem requiring ip to 
be specified at A = 0. To initialize the starting sur- 
face, we need only supply an origin and radius. Taking 
into account that there may be more than one marginally 
trapped surface, it is best to start with an initial surface 
larger than the expected horizon. The values of g-y and 
Kij are required everywhere on the surface to evaluate 
k. These functions can be known analytically or gener- 
ated by evolution codes. As the flow velocity approaches 
zero, F(k — c) — > 0, K ~ c„ and a surface S(c n ) is found 
within a tolerance (e K ). When k = 0, the located surface 
describes a marginally trapped surface. 
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V (K,X) 



nent of the gradient of k with respect to the normals of 
each surface in the level set. The gradient is defined as: 




FIG. 3. Schematic of i/>(c n ), the level set of surfaces in 
We solve for a single surface, S(c n ), in tp(c n ). Multiple levels 
are used in detecting the existence of multiple horizons. 




FIG. 4. Plot of five levels (k = 0.14,0.12,0.10,0.08,0.06) of 
constant divergence of outgoing null geodesies. 

Fig.(|J) shows the level set found in a spacetime 
containing two black holes with coordinate locations 
(-0.954, 0,-0. 3)M and (0.954, 0, 0.3) M. Each 2-surface 
has a constant value of k. We monitor the topology of 
the deforming surface by computing the radial compo- 



ll/tro-i - K n \\ 



(21) 



where ip is the function given in eqn.(|fj). A sharp in- 
crease in the gradient indicates the existence of multiple 
surfaces. To ensure that we do not erroneously abandon 
a single surface, we also monitor the maximum of the Z 2 - 
norm of k. If n is no longer decreasing, we are no longer 
finding a solution to eqn. (Q) ; otherwise the single surface 
is retained. The level flow method is essentially a special 
set of surfaces with properties that let us determine when 
to break. If we only flowed to k = 0, we would not form 
the collection of k ^constant surfaces. 

Once a topology change is indicated, the radii and ori- 
gins for each of the new surfaces are found (note that 
these four parameters for each surface are all that is 
needed to initiate two new trial surfaces). These ori- 
gins and radii are determined using the location of the 
last of the single surfaces. Using this last single surface, 
we can find the origin of the last surface and the location 
on the surface with minimum gradient of the value of k. 
This occurs at the farthest points from the pinching in 
the surface. Picture (a) in fig. (||) shows the last single 
surface with an arrow drawn from the origin to the point 
on the surface with a minimal change in k. The arrow 
indicates a chosen direction. All points lying in this di- 
rection are collected and averaged to find a radius and 
center of mass. All points lying in the opposite direction 
are also collected and used to calculate the radius and 
center of mass for the second surface using the dotted 
arrows in picture (b) of fig.([|). These two sets of radii 
and centers are the initial starting parameters for each 
new trial surfaces. The tracker then flows the two new 
surfaces depicted in (c) of fig.(0) until k = within e K . 




FIG. 5. Two-dimensional schematic representation of the 
three-dimensional decision process to identify the two surfaces 
that will evolve to the two apparent horizons. Fig. (c) depicts 
the two surfaces that will act as new test surfaces. 

The level flow code is only started by the user once, 
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the subsequent flowing to multiple apparent horizons is 
done automatically. 

The advantage to the level flow method is its capabil- 
ity to detect apparent horizons in generic, multiple black 
hole spacetimes from a single reasonable initial guess. 
The drawbacks are the dependence of AA on the spatial 
grid size, AA ~ iV~ 2 where iV 2 = NgN^ is the num- 
ber of grid points, and the fact that we flow to a speed 
of zero (the flow speed approaches zero as k approaches 
zero) . When using apparent horizon tracking in our evo- 
lution code, we will not require knowledge of the apparent 
horizon location to high precision; in fact we can find a 
surface with k < to remove the singularity thus alle- 
viating some of the speed issues. Nonetheless, we plan 
to improve the speed of this algorithm. Some improve- 
ments have already been made to increase the efficiency 
of the current algorithm. The addition of the arctan 2 
function, eqn.(pO|), speeds up the algorithm during the fis- 
sioning process. Monitoring the number of steps needed 
to complete a Crank-Nicholson iteration (see §|v|)) has 
also proven useful in increasing efficiency. 



V. NUMERICAL METHOD 

The previous section described how the level flow 
method is used to solve the apparent horizon equation. 
The resulting parabolic equation is updated using an it- 
erative Crank-Nicholson method updating the variables 
at every A-step. Iterative Crank-Nicholson converges to 
an exact solution of the implicit problem. However, the 
detailed behavior of this convergence jl9| shows that the 
Crank-Nicholson solution at a particular iteration has an 
amplification factor |^4.W| that oscillates around unity. 
The behavior varies in pairs: |.A (n) | < 1 for n = 1,6: 
\A (n) \ > 1 for n = 4,5, etc. while | \A (n) \ - 1| ->■ 
monotonically as n — > oo. n is counting the number of 
iterations it takes to get k = n within the specified tol- 
erance. For the data presented here, a Crank-Nicholson 
iteration of n — 2 or n — 3 was maintained for errors less 
than the spatial grid spacing squared, h 2 . 

The spatial derivatives are approximated to second or- 
der in truncation error using centered finite differencing 
molecules. To verify the convergence of the level flow 
code, we include a plot of the convergence factor versus 
the number of iterations taken in fig. m). 



FIG. 6. Convergence factor for radial variable k with 
h = 0.05, 2h = 0.10, 4ft = 0.20, and AA = 0.0012 for 
Schwarzschild data. Second order accuracy is obtained. 



The convergence factor is given by 



Cf 



&2h — K 4h 
kh - Klh ' 



(22) 



where k is the discretized k and h is the spatial grid spac- 
ing. For a second order scheme, the convergence factor 
in eqn. @ is C f = 4 + 0{h 2 ). 

For the closed form solutions detailed in the next sec- 
tion, the data are given by evaluating the closed form 
analytically on the two-surface. However, the goal is to 
use the level flow method during an evolution including 
evolutions involving a region excised from computational 
consideration. The approach we take to evaluate k from 
a Cartesian grid of data (qg , Kij ) is the same as that used 
and described in Huq [[13| . This approach discretizes the 
apparent horizon equation using Cartesian coordinates 
on 3d-stencils centered on points on the surface. These 
stencils are not aligned with the 3d-lattice from which we 
obtain gij and Kij data. Our apparent horizon surfaces 
are embedded in such lattices and as a result interpola- 
tions must be carried out to obtain the metric data on the 
surface as it evolves. The algorithms and methodology 
for evaluating k are described in detail in [j7 13 1. 



VI. TESTING THE METHOD WITH CLOSED 
FORM SOLUTIONS 



The level flow method of tracking apparent horizons 
has been designed to locate apparent horizons in single 
and multiple black hole spacetimes. To test the level flow 
tracker, we locate apparent horizons in Schwarzschild, 
Kerr, and Brill-Lindquist data. In particular, we also 
demonstrate the level flow method's ability to detect bi- 
nary black holes in the Brill-Lindquist data. 



G 



A. Schwarzschild Data 



B. Kerr Data 



The Kerr-Schild metric provides a closed-form descrip- 
tion of both the Schwarzschild and the Kerr solutions to 
the Einstein equation and is given by: 



9ab ~ Vab + 2Hl a lb, 



(23) 



where r\ a \, is the Minkowski metric, r/ a b =diag(— 1, 1, 1, 1). 
H is a scalar function of the coordinates and l a is an 
ingoing null vector with respect to both the Minkowski 
and full metrics; that is l a satisfies the relation: 



v ab Ub = g ab Ub = 0. 



(24) 



For the ingoing Eddington-Finkelstein form of the 
Schwarzschild solution, the metric given in eqn.(p3|) has 
the scalar function, H, given by: 



H 



M 
r 



and the components of the null vector: 



r 



L= Z - 



(25) 

(26) 
(27) 

(28) 

(29) 



where we have adopted rectangular coordinates (t, x, y, z) 
with r = x 2 + y 2 + z 2 , and M the mass of the black 
hole. 

We track the apparent horizon in this situation for a 
single black hole of mass M. The area and radius of 
the event horizon for the Schwarzschild solution of the 
Kerr-Schild metric is known in closed form |20| given by: 



A = Akt\ 



(30) 



where r + is the event horizon radius and equals 2M . In 
the slicing we have chosen, the apparent horizon coincides 
with the event horizon. Using the level flow method we 
found the apparent horizon to converge to the closed form 
solution giving a 0.35% relative error at a course resolu- 
tion (17 x 17 grid). The area of the tracked apparent 
horizon is computed by 



A r 



VJidxdy, 



(31) 



and converges to the closed form solution, eqn.(|30|). In 
eqn. (|3l|) h is the determinant of the 2-mctric h a i> on the 
apparent horizon surface, and x and y are surface coor- 
dinates. The numerical area is determined from eqn.(|3~l|) 
by calculating the determinant at every point in the grid 
and using a trapezoidal integration scheme p3| . 



The Kerr solution is a second solution given by the 
Kerr-Schild metric, eqn.(|23|). The Kerr solution is the 
solution for a spinning black hole, i.e. a black hole with 
an internal angular momentum per unit mass given by a. 
In rectangular coordinates (t, x, y, z), the scalar function 
and null vector are given by: 

Mr 3 
H = — — — 



and 



(1 



rx + ay ry — ax z 



a* 



a* 



(32) 



(33) 



where /i = (t, x, y, z), M is the mass of the black hole, 
a = J/M is the angular momentum per unit mass of the 
black hole in the z-direction, and r is obtained from: 



x 2 + y 2 



= 1 



(p 2 -a 2 ) + dhp 2 -a 2 ) 2 + a 2 z 2 



(34) 



(35) 



with p = \J x 2 + y 2 + z 2 . The difference here is the addi- 
tion of angular momentum. We test two cases, a = 0.5M 
and a — 0.9M. Fig. (Q) presents a Schwarzschild (a = 
0M) case, together with the a = 0.5M and a = 0.9M 
cases. The solid line is the 6 = tt/2 slice and the dashed 
line is the <p = tt slice. We find the expected result, that 
the deformation in the 6 — n slice increases with a. 




a=0.5M 




FIG. 7. The three plots correspond to the location of the 
apparent horizons for black holes with three different values 
of angular momentum. The units of the graph are M, the 
solid line is the = n/2 slice and the dashed line is the (f> — n 
slice. 
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The radius of the event horizon is given by 



r+ = M + VM 2 - a 2 . 



(36) 



The solution to eqn.(|36j) for a = 0.5M is r + = 1.87M and 
the numerical solution we obtain for the horizon radius 
is r num — 1.87M, with an error of 0.06%. In the a = 
0.9M case, r + = 1.44M and r num = 1.46M, with a 
1.39% error. The area of the horizon for each case can 
be calculated using 



A = 47r(r+ 



(37) 



(generalizing eqn. 



and compared to a numerical 
fl). In the a = 0.5M 



eqn. (pip, computation using eqn. 

case, eqn. (|30|) gives A — 46.89Af z ^numerically we ob- 
tain A num = 46.88M 2 , resulting in a 0.21% error. In 
the a = 0.9M case, A = 36.09M 2 and A num = 36.39M 2 
with a 0.83% error. The errors will decrease as k is driven 
closer to zero. 



a critical separation 1.53(5)M. The apparent horizon at 
the critical separation of 1.535M is shown in fig.(||) using 
the level flow code with 33 2 grid points. 




C. Brill-Lindquist Data 

In this section, we study a binary black hole system 
using Brill-Lindquist data pi] ]. These data are useful 
to us for two reasons: We can verify previous results of 
the critical separation, and study an example of how the 
tracker works in finding multiple apparent horizons. The 
3- metric is time symmetric, K a b — 0, and is conformally 
flat: 



gab = 4> Vab 



where 



= 1 



N 

E 



2n 



(38) 



(39) 



FIG. 8. Separation of 1.535M with a 33^ grid. The area 
was determined numerically to be 184.16M 2 

The horizon found for a separation of d = 1.5M which 
is less than the critical separation, is shown in fig.(H). 
Fig. (|l0|) is a plot of the ^-norm of the maximum of k(0) 
for the separation d = 1.5M at each iteration plotted 
versus the number of A-steps. This is one of the checks 
in the level flow code to ensure that the apparent horizon 
equation is still being solved. We expect the expansion 
to continue to decrease if we have started outside the 
apparent horizon and are flowing inward. As we will see 
in fig.(|l4|), the expansion increases as fission occurs in a 
data set with separated holes. 



and is the number of holes (here N = 2), Mi is the 
mass of the ith black hole and the rj are the radial dis- 
tances from the centers of the black holes 

We use isotropic coordinates to express the metric as 



ds 2 = 4> 4 (dr 2 +r 2 d9 2 



2 ■ 2 

r sm 



with 



2d,rcos( 



(40) 



(41) 



where di is the distance between the holes and the center 
of the coordinate system. When they are far apart, each 
hole has an apparent horizon of radius M/2 in these co- 
ordinates. The area of each of the holes when they are 
well separated is 16irM 2 . 

The limiting separation of the holes between single and 
double horizons was found by Brill and Lindquist ]2l| ] to 
be 1.56M, Cadez 1.534M±0.002Af ||, [1.5M, 1.6M] by 
Alcubierre et al. ||, and 1.535M by Huq ||. We found 



FIG. 9. Separation of 1.5M with a 33^ 
185. 41M 2 . 



grid. The area is 
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FIG. 11. The starting surface of the level flow method for 
a separation of d=2.0M 



FIG. 10. The absolute value of the maximum of expansion, 
k, per iteration, A every 20th step. The kinks at A = 1000 
and A = 2000 are from restarting the code with a different 
A-step. 

As we increase the separation between the two holes to 
a separation greater than the critical separation, we can 
test the apparent horizon tracker in the case of multiple 
apparent horizons. We demonstrate with a separation 
of d = 2.0M. The initial surface flows to the point of 
fissioning where the topology of the surface changes from 
a one surface into two surfaces. Fig. (11) is a plot of 
the initial surface that begins the flow. The level set 
found during this flow is depicted in fig. (|l^) . Each of the 
surfaces in fig.(|l2|) has a constant expansion, k = c n and 
was used to indicate a topology change in the test surface. 
The values for the expansion are c\ — 0.14, C2 = 0.12, 
c 3 = 0.1, and C4 = 0.08. The last single surface just 
before the topology change is not a surface in the level 
set; it is plotted in fig. (|13|) . At this point the tracker 
begins to flow two surfaces. 




FIG. 12. The level set of surface for the d=2.0M case. Each 
surface has a constant re = c„ at each point. In this case 
the values for c„ are ci = 0.14, C2 = 0.12, C3 = 0.1, and 
C4 = 0.08. The level set is used to indicate the change in 
topology associated with multiple surfaces. 

In contrast to a separation of d = 1.5M where there 
is no fission, here as fissioning becomes imminent, the k 
begins to increase. Fig.(|l4|) is a plot of the absolute value 
of the maximum across the surface of the expansion, k, 
versus A up to the point of fission. The increase in the 
expansion is one of the signals of imminent fission. As the 
algorithm tries to find a surface with n = everywhere, 
it is driven into two surfaces. Once the new surfaces are 
found, the maximum of the expansion begins a monotonic 
decrease as in fig. (10). 




FIG. 13. The single surface is about to fission into two 
surfaces. 
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d = 3.0M 




50 100 150 

X 



FIG. 14. The absolute value of the maximum of expansion, 
K, per iteration is plotted. The increase in the expansion is 
caused by imminent fission. 

The exaggerated peanut shape in fig.jHf) and fig.fl6|) 
is taken for the same A- value as fig . (p~3|) . 

Once the fissioning is detected by the code, it auto- 
matically begins flowing two new surfaces of the same 
resolution as the parent surface. The series of snapshots 
shown in fig.(^5|) and fig. (|l6|) is a subset of the set of 
surfaces found by the apparent horizon tracker as it fol- 
lows the fission of the trial surface into two surfaces. The 
tracker starts with a spherical starting surface that de- 
forms along the gradient field. 




FIG. 15. This series of snapshots depicts the flow of an 
initial surface until its fission for the binary Brill-Lindquist 
black holes separated by 2M. The lower left plot is a first 
try at determining the final two surfaces. The cusps are due 
to a typical drawback associated with using points to define 
the flowing surface. The points crowd together in regions of 
greater flow. The next snapshot, on the lower right, shows the 
code's automatic correction; and shows the apparent horizons 
of the binary Brill-Lindquist data to an accuracy of 10 -4 . 

As the points defining the surface flow, the distance 
between the points can become too small for the finite 
difference scheme at that resolution. Redistribution of 
the points on the surface is taken care of automatically 
by updating the center and radius. 




-2-1012 
x[M] 



FIG. 16. The series of pictures shown in fig.(|l5|) are placed 
in one plot. The outer surface is the initial guess, the "peanut 
surface" is the surface that is found indicating the need to 
search for two surfaces, and the inner approximate spheres 
result from locating the apparent horizons for Brill-Lindquist 
data. 



VII. APPARENT HORIZONS IN A GRAZING 
COLLISION 

As stated in the Introduction, one of the main moti- 
vators of this work is to have an apparent horizon finder 
that can locate disjoint horizons during the evolution. 
This entails 1) finding the horizons without a good ini- 
tial guess, and 2) detecting the topology change from two 
disjoint horizons to one horizon. To demonstrate the level 
flow's ability to carry out 1) and 2), we report the results 
of apparent horizon tracking in the particular case of a 
short evolution of two spinning, Kerr-Schild black holes 
using the Binary Black Hole Grand Challenge Alliance 
Cauchy code [24]] . A future paper ]25|] will discuss the 
details of the evolution. 
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The evolution is free, i.e. the momentum and Hamil- 
tonian constraint equations are only used as checks dur- 
ing the evolution. Since we cannot hold infinity on the 
grid in this formalism, we must specify outer boundary 
conditions for the dynamic variables, g^ and Kij. For 
the following work, we specify analytic outer boundary 
conditions with blending |26| between the analytic and 
numeric regions. 

To specify initial data for two spinning, boosted black 
holes we use superposed Kerr-Schild black holes. We 
chose a Kerr-Schild metric |20| for two reasons: 1) The 
metric is well defined at the event horizon, and 2) The 
metric is Lorentz form-invariant in a simple sense, under 
boosts (v > 0). The superposed data are constructed 
in an approximate manner by a conformal method based 
on the superposition of two isolated, boosted Kerr black 
holes. 

The initial data follows from Matzner et al. |27| and 
was first implemented by Correll The data is the 

superposition of two, isolated, boosted Kerr-Black holes 
with individual metrics given by eqn. (p3|) . The resulting 
superposed metric is: 

9ij = {i)9H + (2)9ij ~ mj (42) 

with the symbol indicating a quantity conformally 
related to the physical metric, gij = <fi gij. 

(1) 9ij = 'nij + (i)H{ri)(i)k(i)lj and (43) 

(2) 9ij = Vij + (2)H(r 2 )( 2 )h(2)lj (44) 

are the the isolated Kerr-Schild metric forms with /j and 
H corresponding to the single black holes. The two holes 
have comparable masses, Mi ~ M 2 , coordinate separa- 
tion ri2, and velocities V\ and v 2 assigned to them. For 
the argument of H and lj , we use 

ri 2 = (x - XxYix - XiySij and (45) 
r 2 2 = (x-x 2 y{x-x 2 yS t] (46) 

with x\ and x 2 J the coordinate positions of the holes on 
the initial slice. 

The extrinsic curvatures of the two isolated black holes 
are added to obtain a trial K a b for the binary black hole 
system given as: 

(Q)Kij = (i)Kij + ( 2 )Kij. (47) 

The subscript indicates that this is an approximation 
to the true extrinsic curvature of the binary black hole 
spacetime. (i)Kij and ( 2 )Kij are the individual extrin- 
sic curvatures associated with the isolated Kerr-Schild 
metric and their indices are raised and lowered by their 
individual metrics, eqn. ([43]). 

For the data we describe here (holes center initially 
separated by a coordinate distance exceeding 10M where 
M is the mass of one of the holes), we expect that an 
initial value solution will be ps 10% in error on the domain 



outside of the excision volume. See further discussion in 
p{|. We set the lapse function to: 

a = ai + a 2 — 1, 
and the shift vector to: 

The run presented in this paper has a grid 81 x 81 x 
81 in Cartesian coordinates (x, y, z) with a domain of 
(±10M, ±10M, ±10M) resulting in a spatial resolution 
of M/4. The data represent two black holes in a grazing 
collision. The holes are set initially at (5M, IM,0M) 
and (— 5Af, — 1M, 0M) in Cartesian coordinates with a 
boost speed of ±0.5cc toward one another and each has 
an angular momentum per unit mass of a = 0.5M in 
the (— )z-direction. Fig.(|l7j) is the initial configuration of 
this run; note that a naive sum of the spin and the orbit 
angular momentum yields zero for this configuration. 

We post-process the data obtained from the evolu- 
tion. For the purposes of this paper, we track the appar- 
ent horizons at three specific times during the evolution, 
namely t = OA/, 2.8A/, and 3.4M. At 0M, the apparent 
horizons of the initial data are found, at 2.8M two dis- 
joint apparent horizons are found; and finally, at 3.4M a 
single apparent horizon is found. For the horizons shown 
here, the level flow method used a sphere of radius 8M to 
initialize each run. Fig. [l8] is a plot of the horizons with 
time going up the page. The lowest plot is of t = 0M, 
with each horizon being a sphere centered at coordinates 
(±5M, ±1A/,0A/). The middle plot shows the horizons 
at a later time, t = 2.8M. Here the deviation in shape as 
the horizons accelerate towards each other is seen. The 
final plot at the top of fig. [l8] is the first single apparent 
horizon that envelops both black holes at t = 3.4M . 

The areas for the apparent horizons at t = OAf are 
A = 43.6A/ 2 for each hole. At t = 2.8M, the horizons 
have deviated from a spherical shape and the areas for 
each hole are A = 44.2M 2 , giving a measure of the ac- 
curacy to which we can maintain their areas constant. 
The area of the merged apparent horizon at t — 3.4M is 
Amerged = 184M 2 . According to the black hole area the- 
orem of Hawking and Ellis pp[ , the area of the merged 
event horizon must equal at least the combined area of 
the individual event horizons. Although no strong state- 
ments can be made about the area of an apparent hori- 
zon, we do find A merge( i > A\ + A 2 in a consistent man- 
ner. We can further surmise that the final maximum area 
we could expect based on the initial configuration should 
be that of a Schwarzschild black hole of mass 2Af , giving 
an area of approximately 201Af 2 . In some sense the area 
predicted by the Schwarzschild case is an upper bound. 
We see a 8.5% deviation from that "idealized" case. In 
view of this upper limit, 8.5% may be an indication of 
the greatest amount of gravitational radiation up to the 
time of merger (t = 3.4A/) given our approximate initial 
data, gauge condition, and boundaries. 
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VIII. CONCLUSION 



a=-0.5z 



y[M] 



v=0.5x 



x[M] 



FIG. 17. The configuration of the initial data for the graz- 
ing collision. The initial mask position is indicated by the 
"circle" centered on each hole. The angular momentum per 
unit mass and velocity of each of the holes is also represented. 
The total angular momentum (spin plus orbit) of the initial 
configuration is zero. 




Apparent horizon location and tracking constitute an 
important part of numerical evolutions of black hole 
spacetimes using excision techniques. We have demon- 
strated a method for finding apparent horizons in situa- 
tions where the location of the apparent horizon may not 
be known; hence a good initial guess for the finder may 
not be possible. The method we have discussed works 
with generic 3-metric and extrinsic curvature data and 
with an arbitrary initial starting surface. Furthermore, 
the method is capable of detecting a topology change 
as the finder flows towards the apparent horizon. This 
ability is important for situations where there are mul- 
tiple apparent horizons in the data. This allows us to 
locate apparent horizons in binary black hole evolutions 
without knowing where the apparent horizons are; and it 
allows us to locate the first single apparent horizon that 
forms at the merger of two black holes. The level flow 
method is successful at locating the apparent horizons in 
generic spacetimes as demonstrated by the Schwarzschild 
and Kerr data. It also found multiple apparent hori- 
zons starting from a single starting surface as demon- 
strated with the Brill-Lindquist data. Most importantly, 
the level flow method has been successful at identifying 
apparent horizons in a binary black hole evolution involv- 
ing two Kerr-Schild black holes. Beginning with a single 
guess surface, two discrete apparent horizons were found 
at early times, and the later single merged horizon was 
found. 

One of the drawbacks of this method currently is its 
slow convergence property due to the parabolic nature 
of the equation solved. This, however does not pose a 
problem since the level-flow method can be used in con- 
junction with other methods which may be more efficient 
given a good initial guess. The level-flow method has the 
definite advantage of being capable of finding multiple 
surfaces in the data. It can be used to get extremely 
good initial guesses for other methods that converge more 
quickly close to the solution. We are currently using the 
level flow method in this manner in numerical evolutions 
of black hole collisions. 




FIG. 18. The apparent horizons are plotted with the evo- 
lution time increasing up the page. The times plotted are 
0.0M, 2.8M, and 3AM. 
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Generic Tracking of Multiple Apparent Horizons with Level Flow 
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We report the development of the first apparent horizon locator capable of finding multiple apparent 
horizons in a "generic" numerical black hole spacetime. We use a level-flow method which, starting 
from a single arbitrary initial trial surface, can undergo topology changes as it flows towards disjoint 
apparent horizons if they are present. The level flow method has two advantages: 1) The solution 
is independent of changes in the initial guess and 2) The solution can have multiple components. 
We illustrate our method of locating apparent horizons by tracking horizon components in a short 
Kerr-Schild binary black hole grazing collision. 



I. INTRODUCTION 

Our goal is to investigate the strong field regime of gen- 
eral relativity. In particular we wish to focus on the study 
of coalescing black hole binaries. Over the last three 
decades since the pioneering work of Cadez, Smarr, Ep- 
pley and others, advances in computing technology, nu- 
merical algorithms and techniques and our understand- 
ing of the underlying physics have advanced to a point 
where we are able to carry out simulations of binary black 
hole collisions in 3+1 dimensions. One of the outcomes 
of such simulations will be an understanding of the un- 
derlying physics of the problem; and, therefore, a pre- 
diction and understanding of the gravitational radiation 
content. A detailed knowledge of how the resultant grav- 
itational waveforms relate to the physical parameters of 
the binaries that produce them will be of importance to 
gravitational wave observatories (such as LIGO, VIRGO, 
TAMA, GEO600) now under construction around the 
world. With the building of these observatories, we stand 
at the epoch of the first direct observations of astrophys- 
ical sources that involve strong field general relativity. 

The orbit and merger of two black holes is one candi- 
date source for ground based detection of gravitational 
waveforms. This is of great interest to the relativity com- 
munity. The binary black hole problem is a two-body 
problem in general relativity. It is a stringent dynam- 
ical test of the theory. However, studying spacetimes 
containing multiple black holes involves solving the Ein- 
stein equations, a complex system of non-linear, dynamic, 
elliptic-hyperbolic equations intractable in closed form. 
The intractability of the problem has led to the develop- 
ment of numerical codes capable of solving the Einstein 
equations. 

Our approach to numerically solving the Einstein equa- 
tions involves reformulating them as an initial value prob- 
lem. In this 3+1 formulation [1], spacelike hypersurfaces 
parameterized in time foliate the spacetime. The result- 
ing equations are coupled elliptic and hyperbolic differ- 
ential equations of the 3-metric, gij, and the extrinsic 



curvature, Kij. The initial value problem is solved by 
specifying a hypersurface at an instant of time, say t = 0, 
and evolving to the next hypersurface at time t = St with 
the evolution equations to obtain g^ and at the next 
time t = St. 

One vital issue in numerically solving the Einstein 
equation describing spacetimes containing black holes is 
handling the physical singularities. As one approaches 
the singularity, the values of the fields being computed 
approach infinity; therefore, a region containing the sin- 
gularity must be avoided to keep the computation from 
halting. Excision techniques are promising in avoiding 
the singularity. Excising the singularity involves locat- 
ing a region interior to the event horizon containing the 
singularity on each evolving hypersurface. This region is 
then "masked" from the computation. The derivatives at 
the boundary between the masked region and the compu- 
tational domain are handled using causal differencing, a 
differencing scheme [2] that respects the causal structure 
of the spacetime. 

In deciding where in space to excise we use the appar- 
ent horizon as opposed to the event horizon. By its very 
nature, the event horizon is a global construct depending 
on the entire spacetime. The apparent horizon, a local, 
i.e. spacelike 2-surface is more suitable. Following a sug- 
gestion by Unruh [3], the apparent horizon is used to 
define the excised region to be masked at each time dur- 
ing the evolution. Apparent horizons are defined locally 
for each time, and exist at, or inside of, the event horizon. 
In some spacetimes, choices of foliation may lead to the 
absence of an apparent horizon. When the discussions 
in this refer to a hypersurface, it is assumed an apparent 
horizon exists on that hypersurface. 

Recent three-dimensional horizon locator codes are ca- 
pable of finding the location of an apparent horizon 
in generic single black hole spacetimes [4-9] and two 
[10,11] are capable of finding disjoint multiple apparent 
horizons in the special case of conformally flat binary, 
time-symmetric black hole spacetimes. Multiple appar- 
ent horizon finding algorithms will be necessary in sim- 
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ulations of generic binary black hole spacetimes. The 
method presented in this paper, called the level flow 
method, is capable of detecting multiple apparent hori- 
zons in generic spacetimes. The level flow algorithm has 
two advantages: 1) It is independent of a good initial 
guess and 2) It is capable of following the surface through 
a change in topology. In level flow, the apparent horizon 
equation is reformulated as a parabolic equation and a 
set of surfaces are flowed with speeds dependent on the 
expansion of the outgoing null vectors normal to each 
surface. 

The purpose behind developing the level flow method 
of tracking apparent horizons is to have a method capa- 
ble of detecting multiple apparent horizons on any given 
hypersurface without a good guess. Specifically, we want 
a tracker that can detect the transition from a double to 
a single apparent horizon in single time step without a 
prior knowledge of the transition. 

In the rest of this paper we discuss apparent horizons 
in general and current 3D work in § 77 and § III. In 
§ IV we describe the level flow method in detail and give 
a brief description of the numerical method involved in 
solving the apparent horizon is in § V. We demonstrate 
the use of the level flow tracker on model data and a 
binary black hole grazing collision in § VI and § VII. 

II. APPARENT HORIZONS 




The expansion of the outgoing normals, V fc°, can be 
rewritten in terms of quantities defined on the hypersur- 
face: 

k = D a s a -K + s a s b K ab . (3) 

D a denotes covariant differentiation with respect to g a b, 
K a i, is the extrinsic curvature, and K is g ab K a i,. In fact, 
there is a level set of surfaces in £ parameterized by k. 
Each surface in the level set is defined by the constant 
expansion of its null vectors such that 

k = c„, (4) 

where c„ are constants labeled by the positive integer n. 
Marginally trapped surfaces are members of this set for 

k = 0. (5) 

Eqn.(5) is called the apparent horizon equation since the 
apparent horizon is the outermost surface with k = in 
S. 

The S 2 topology of the apparent horizon naturally 
lends itself to characterization via spherical coordinates. 
The function, 

iP = r-h(6,cf>) (6) 

is a level set of 2-spheres in S, where h(6,<j>) is called the 
apparent horizon shape function. A marginally trapped 
surface has ip = 0. The spacelike normals to S are defined 
from eqn. (6) such that 

s* = gVdrf/y/jjudidrfW (7) 

is the spacelike vector field at every point of S. 
The apparent horizon equation in spherical coordinates 
(h(8, 4>),0, 4>) is a 2-dimensional problem in 8 and <j>. 



FIG. 1. Representation of a 2-sphere embedded in a hyper- 
surface, E. 

M is the spacetime with metric A g a b foliated by hyper- 
surfaces £ parameterized by t with 3-metric g a b. Let S 
be a surface with S 2 topology on £. The apparent hori- 
zon is the outermost marginally trapped surface in S, a 
surface with zero expansion. The zero expansion of the 
surface is defined in terms of outgoing null vectors to S, 
k a , such that k a have zero divergence 

V a fc a = 0, (1) 

where V a is the covariant derivative associated with 4 g a b- 
Referring to fig. (1), k a is defined in terms of the space- 
like normals to S, s a , and the future directed timelike 
normals, n a , such that: 

k a = s a + n a . (2) 



III. CURRENT 3-D METHODS 

The approaches to solving the apparent horizon equa- 
tion on a three-dimensional hypersurface can be ad- 
dressed roughly in two categories: methods that solve 
the apparent horizon equation directly and methods that 
solve it by first recasting it as a parabolic equation. This 
paper does not address spherical and axi-symmetric ap- 
proaches. One of the first three-dimensional apparent 
horizon solvers was published by Nakamura, Kojima, and 
Oohara [9]. They directly solve the apparent horizon 
equation by using spherical harmonic basis functions to 
expand the apparent horizon shape function, h(9, </>) into 

h{6,<t>) = £ Y, a imYi m (0,<f>). (8) 

/=0 m=-l 

This method is called the pseudospectral method. A fi- 
nite number of the coefficients, {a; TO } parameterize the 
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horizon shape function, and the maximum l max depends 
on the computation and deviations from a sphere. The 
apparent horizon equation can then be solved by writing 
it as 



\n{ai r 



0, 



(9) 



and using functional integration routines to find the co- 
efficients ai m . Others have used similar methods [5,8,4]. 

In another approach to direct solutions of the appar- 
ent horizon equation is to treat it as a boundary value 
problem. One notes that a discretization of this equa- 
tion leads to a system of algebraic equations which can 
then be solved via Newton's method. Thornburg [12] dis- 
cusses applications of Newton's method to this problem 
in general and shows results in axisymmetry. Huq [13] 
has implemented a similar algorithm based on Newton's 
method that utilizes Cartesian coordinates to difference 
the equations. 

Tod [14] first suggested the use of curvature flow in 
solving the apparent horizon equation by recasting it as 
a parabolic equation. Bernstein [15] implemented Tod's 
suggestion in axisymmetry. Gundlach [6] introduced a 
fast flow method which combines the ideas of the flow 
method with the pseudospectral method. Pasch [11] and 
Diener [10] implemented a similar method, a level-set 
method, in three-dimensions and found discrete apparent 
horizons in multi-black-hole spacetimes; however these 
spacetimes were confined to be conformally flat and time- 
symmetric. 

Each of the approaches briefly described above, solving 
the apparent horizon equation directly or solving it via 
a parabolic equation, has its advantages. Direct meth- 
ods tend to be faster while flow methods do not rely on 
"good" initial guesses. However, none of these methods 
are applicable to the generic, multi-black-hole problem. 
Herein lies the motivation behind the level flow method. 
The level flow method is the only method designed to 
locate discrete apparent horizons in generic spacetimes 
containing one or two discrete horizons. 



IV. LEVEL FLOW METHOD 

A. Curvature Flow 

The level flow method is a hybrid flow/level— set 
method. The previous section mentioned the flow 
method, this section gives more detail on the flow method 
which is the foundation of the level flow method. The 
flow approach, as suggested by Tod, is to rewrite the 
apparent horizon equation as the speed function in a 
parabolic equation. In the case of a time symmetric 
hypersurface, in which K a b = 0, the apparent horizon 
equation reduces to the condition for a minimal surface, 
D a s a = 0. In this case, the surface 5 is at a local ex- 
tremum of the area. The starting surface, S(X = 0), is 



parameterized by coordinates x a and evolved in terms of 
a parameter A. The equation of motion is 



dx a 
~8X 



-Hs a 



(10) 



where dx a /d\ is a vector field, and H is the mean curva- 
ture, which is the trace of extrinsic curvature associated 
with embedding S in £ given by 



H = D a s a . 



(11) 



Eqn.(10) is the gradient flow for the area functional. The 
area decreases monotonically with increasing A. Grayson 
[16] has shown that a surface deforming under its gra- 
dient field (Eqn.(10)) will evolve to a stable minimum 
surface (surface is local minimum of area) if there is one, 
otherwise to a point. 

In numerical relativity, we are interested in the generic 
case, with K a b ^ 0, for which the marginally trapped 
surfaces differ from minimal surfaces, the surfaces are not 
extrema of the area. However, Tod suggests an equation 
similar to Eqn.(10) as a curvature flow: 



dx a 
~8X 



F( K )s a 



(12) 



using F(k) = —k where k = D a s a + s a s h K a i > — K as 
in eqn.(3). We have found eqn. (12) to be a successful 
practical implementation of the flow method. 



B. Level Flow 

Eqn. (12) gives us an initial value problem. Given in- 
formation about the system at some initial A, eqn. (12) 
will describe the system for all its future propagation in 
A. Directly solving eqn. (12) will lead to a successful de- 
tection of single apparent horizons; however, solving it 
directly does not ensure correct handling of a topology 
change which is necessary for detection of multiple appar- 
ent horizons. By combining the flow method with a level- 
set idea however, this topology change can be effected 
and multiple apparent horizons can be tracked starting 
from a single guess surface. 

First eqn. (12) is recast from an equation governing 
the motion of the coordinates parameterizing S, namely 
x a , to an equation governing the motion of the surface 
ip. Noting that ip is parameterized by A, 



dtp 
dX 



dx a dip 
dX dx a 



(13) 



by the chain rule, and multiplying eqn. (12) by J^L gives 
dtb 



8X 



F( K )s a dip/dx a . 



(14) 



Using 
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s a =9 ab -^ b IW\\ (15) 

and 

11^11= i/ii^. el 

the test surface's flow is given by: 

^=F( K )||W||. (17) 

Eqn.(17) is a reformulation of eqn. (12) and will flow the 
surface, ip, to a marginally trapped surface at tp = when 

K = 0. 

The strength of the flow methods is their ability to lo- 
cate a surface with k = given any non-pathological 
initial surface. For example, the apparent horizon in 
a spherically symmetric spacetime (Schwarzschild) was 
found by flowing an initial surface shaped as a leaf, see 
Fig. (2). This ability is especially important when track- 
ing horizons during evolutions of binary black hole space- 
times. In this case, finding apparent horizons for two dis- 
crete apparent horizons in each £(i) involves flowing two 
initial guesses, one for each horizon [17]. On £(£ = 0), the 
location of the apparent horizons may be known; how- 
ever, as the black holes accelerate the task of guessing 
the locations of the two horizons becomes more difficult. 
Further, the two horizons merge into a single horizon at a 
single instant of time, rendering a good initial guess diffi- 
cult. Some way of determining when two horizons merge 
into a single horizon is necessary. The level flow method 
takes care of these issues by not requiring a good ini- 
tial guess (ip(\ = 0)) and by detecting multiple apparent 
horizons from a single guess (ip(\ = 0)). 
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FIG. 2. Example of a surface undergoing flow on a spher- 
ically symmetric hypersurface. The figure labeled (a) is 
ip(\ = 0), the initial guess. Figure (b) is ip{\) representing 
the final solution, the apparent horizon. 

The level flow method differs from the flow method in 
the specification of the speed function, F(k). F(k) = 
determines when the propagation of the surface stops. 
The flow method is in the form of eqn. (12), in which 
F(k) = —k. A good choice since F(k) = —k = indi- 
cates a marginally trapped surface; but this choice will 
not flow ip though a fission. In general the scheme fails 
as the surface pinches due to ill-defined normals at the 
surface. The level flow method alleviates this problem by 
looking for indications that the surface topology is about 
to change before the pinching occurs. (Another method 
which was introduced by Sethian and Osher [18] for non- 
relativistic problems is to flow a higher dimensional sur- 
face in which tp is embedded. This higher dimensional 
surface does not fission. This has only been implemented 
in a time-symmetric spacetime [11] and requires more 
computational power due to the extra dimension in the 
problem.) 

The level flow method flows a set of two-dimensional 
surfaces in £ parameterized by k. We call the set of 
surfaces a level set and label the set S(c„). Each surface 
has a constant value of k = c„ everywhere on it. The set 
of surfaces is defined by varying c n as the flow progresses 

c„+i = c„ ± Ac, (18) 

where (+) indicates outward flow, (— ) inward flow, and 
Ac oc 1 1 k ||2. Each surface obeys the equation of motion 
given in eqn. (17) with F(k) defined to flow to multiple 
surfaces. We choose two options for the speed function: 

F(k) = K-c n (19) 
F(k) = (k - Cn) arctan 2 (-^^). (20) 

K a 

As k — c n — > 0, both functions are solving for a partic- 
ular surface in the level set, S(c n ). The second func- 
tion, eqn. (20), behaves similarly to the first but allows 
for larger time steps near a fissioning surface because it 
moves points further from the k — c n = surface faster 
than the points closer to this surface. 

Eqn. (17) is an initial value problem requiring ip to 
be specified at A = 0. To initialize the starting sur- 
face, we need only supply an origin and radius. Taking 
into account that there may be more than one marginally 
trapped surface, it is best to start with an initial surface 
larger than the expected horizon. The values of <?y and 
Kij are required everywhere on the surface to evaluate 
k. These functions can be known analytically or gener- 
ated by evolution codes. As the flow velocity approaches 
zero, F(k — c) — > 0, k ~ c n and a surface S(c n ) is found 
within a tolerance (e K ). When k = 0, the located surface 
describes a marginally trapped surface. 
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nent of the gradient of k with respect to the normals of 
each surface in the level set. The gradient is defined as: 




FIG. 3. Schematic of tp(c n ), the level set of surfaces in £(t). 
We solve for a single surface, S(c n ), in ip(c n )- Multiple levels 
are used in detecting the existence of multiple horizons. 




FIG. 4. Plot of five levels (k = 0.14,0.12,0.10,0.08,0.06) of 
constant divergence of outgoing null geodesies. 

Fig. (4) shows the level set found in a spacetime 
containing two black holes with coordinate locations 
(-0.954,0, -0.3)M and (0.954,0, 0.3)M. Each 2-surface 
has a constant value of k. We monitor the topology of 
the deforming surface by computing the radial compo- 



| K n-l 



tpn\ 



(21) 



where ip is the function given in eqn.(6). A sharp in- 
crease in the gradient indicates the existence of multiple 
surfaces. To ensure that we do not erroneously abandon 
a single surface, we also monitor the maximum of the li- 
norm of k. If k is no longer decreasing, we are no longer 
finding a solution to eqn.(4); otherwise the single surface 
is retained. The level flow method is essentially a special 
set of surfaces with properties that let us determine when 
to break. If we only flowed to k = 0, we would not form 
the collection of k =constant surfaces. 

Once a topology change is indicated, the radii and ori- 
gins for each of the new surfaces are found (note that 
these four parameters for each surface are all that is 
needed to initiate two new trial surfaces). These ori- 
gins and radii are determined using the location of the 
last of the single surfaces. Using this last single surface, 
we can find the origin of the last surface and the location 
on the surface with minimum gradient of the value of k. 
This occurs at the farthest points from the pinching in 
the surface. Picture (a) in fig. (5) shows the last single 
surface with an arrow drawn from the origin to the point 
on the surface with a minimal change in k. The arrow 
indicates a chosen direction. All points lying in this di- 
rection are collected and averaged to find a radius and 
center of mass. All points lying in the opposite direction 
are also collected and used to calculate the radius and 
center of mass for the second surface using the dotted 
arrows in picture (b) of fig. (5). These two sets of radii 
and centers are the initial starting parameters for each 
new trial surfaces. The tracker then flows the two new 
surfaces depicted in (c) of fig. (5) until k = within e K . 




FIG. 5. Two-dimensional schematic representation of the 
three-dimensional decision process to identify the two surfaces 
that will evolve to the two apparent horizons. Fig. (c) depicts 
the two surfaces that will act as new test surfaces. 

The level flow code is only started by the user once, 
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the subsequent flowing to multiple apparent horizons is 
done automatically. 

The advantage to the level flow method is its capabil- 
ity to detect apparent horizons in generic, multiple black 
hole spacetimes from a single reasonable initial guess. 
The drawbacks are the dependence of AA on the spatial 
grid size, AA ~ N~ 2 where N 2 = NgN^ is the num- 
ber of grid points, and the fact that we flow to a speed 
of zero (the flow speed approaches zero as k approaches 
zero) . When using apparent horizon tracking in our evo- 
lution code, we will not require knowledge of the apparent 
horizon location to high precision; in fact we can find a 
surface with k < to remove the singularity thus alle- 
viating some of the speed issues. Nonetheless, we plan 
to improve the speed of this algorithm. Some improve- 
ments have already been made to increase the efficiency 
of the current algorithm. The addition of the arctan 2 
function, eqn.(20), speeds up the algorithm during the fis- 
sioning process. Monitoring the number of steps needed 
to complete a Crank-Nicholson iteration (see §V)) has 
also proven useful in increasing efficiency. 

V. NUMERICAL METHOD 

The previous section described how the level flow 
method is used to solve the apparent horizon equation. 
The resulting parabolic equation is updated using an it- 
erative Crank-Nicholson method updating the variables 
at every A-step. Iterative Crank-Nicholson converges to 
an exact solution of the implicit problem. However, the 
detailed behavior of this convergence [19] shows that the 
Crank-Nicholson solution at a particular iteration has an 
amplification factor \A^\ that oscillates around unity. 
The behavior varies in pairs: |_4(")| < 1 for n = 2,6; 
\A (n) \ > 1 for n = 4,5, etc. while | \A (n) \ — 1 1 — > 
monotonically as n — > oo. n is counting the number of 
iterations it takes to get k = k within the specified tol- 
erance. For the data presented here, a Crank-Nicholson 
iteration of n = 2 or n = 3 was maintained for errors less 
than the spatial grid spacing squared, h 2 . 

The spatial derivatives are approximated to second or- 
der in truncation error using centered finite differencing 
molecules. To verify the convergence of the level flow 
code, we include a plot of the convergence factor versus 
the number of iterations taken in fig. (6). 




FIG. 6. Convergence factor for radial variable n with 
h = 0.05, 2h = 0.10, Ah = 0.20, and AA = 0.0012 for 
Schwarzschild data. Second order accuracy is obtained. 

The convergence factor is given by 

C f = k2h ~ kih , (22) 

Kh - K2h 

where k is the discretized k and h is the spatial grid spac- 
ing. For a second order scheme, the convergence factor 
in eqn. (22) is C f =4 + 0(h 2 ). 

For the closed form solutions detailed in the next sec- 
tion, the data are given by evaluating the closed form 
analytically on the two-surface. However, the goal is to 
use the level flow method during an evolution including 
evolutions involving a region excised from computational 
consideration. The approach we take to evaluate k from 
a Cartesian grid of data (gij,Kij) is the same as that used 
and described in Huq [13]. This approach discretizes the 
apparent horizon equation using Cartesian coordinates 
on 3d-stencils centered on points on the surface. These 
stencils are not aligned with the 3d-lattice from which we 
obtain gij and Kij data. Our apparent horizon surfaces 
are embedded in such lattices and as a result interpola- 
tions must be carried out to obtain the metric data on the 
surface as it evolves. The algorithms and methodology 
for evaluating k are described in detail in [7,13]. 

VI. TESTING THE METHOD WITH CLOSED 
FORM SOLUTIONS 

The level flow method of tracking apparent horizons 
has been designed to locate apparent horizons in single 
and multiple black hole spacetimes. To test the level flow 
tracker, we locate apparent horizons in Schwarzschild, 
Kerr, and Brill-Lindquist data. In particular, we also 
demonstrate the level flow method's ability to detect bi- 
nary black holes in the Brill-Lindquist data. 
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A. Schwarzschild Data 



B. Kerr Data 



The Kerr-Schild metric provides a closed-form descrip- 
tion of both the Schwarzschild and the Kerr solutions to 
the Einstein equation and is given by: 



9ab — f]ab + 2Hl a l b , 



(23) 



where rj a b is the Minkowski metric, rj a b =diag(— 1, 1, 1, 1). 
H is a scalar function of the coordinates and l a is an 
ingoing null vector with respect to both the Minkowski 
and full metrics; that is l a satisfies the relation: 



v ab Ub 



g ab Ub = 0. 



(24) 



For the ingoing Eddington-Finkelstein form of the 
Schwarzschild solution, the metric given in eqn.(23) has 
the scalar function, H, given by: 



H 



M 
r 



and the components of the null vector: 



1 

x 

r 

y 

r 

z 



(25) 

(26) 
(27) 

(28) 

(29) 



where we have adopted rectangular coordinates (t, x, y, z) 
with r = \J x 2 + y 2 + z 2 , and M the mass of the black 
hole. 

We track the apparent horizon in this situation for a 
single black hole of mass M. The area and radius of 
the event horizon for the Schwarzschild solution of the 
Kerr-Schild metric is known in closed form [20] given by: 



Auri 



(30) 



where r + is the event horizon radius and equals 2M. In 
the slicing we have chosen, the apparent horizon coincides 
with the event horizon. Using the level flow method we 
found the apparent horizon to converge to the closed form 
solution giving a 0.35% relative error at a course resolu- 
tion (17 x 17 grid). The area of the tracked apparent 
horizon is computed by 



A r 



/ Vhdxdy, 
Js 



(31) 



and converges to the closed form solution, eqn.(30). In 
eqn.(31) h is the determinant of the 2-metric h a b on the 
apparent horizon surface, and x and y are surface coor- 
dinates. The numerical area is determined from eqn.(31) 
by calculating the determinant at every point in the grid 
and using a trapezoidal integration scheme [13]. 



The Kerr solution is a second solution given by the 
Kerr-Schild metric, eqn.(23). The Kerr solution is the 
solution for a spinning black hole, i.e. a black hole with 
an internal angular momentum per unit mass given by a. 
In rectangular coordinates (t,x,y,z), the scalar function 
and null vector are given by: 

Mr 3 

H = -^-^ (32) 



2 r 2 



and 



(1 



r 2 + a 2 z 



rx + ay ry — ax z, 
r' 2 + a 2 ' r 2 + a 2 ' r ' 



(33) 



where p = (t,x,y,z), M is the mass of the black hole, 
a = J/M is the angular momentum per unit mass of the 
black hole in the z-direction, and r is obtained from: 



■x 2 + y 2 z 2 

— + — 

2 + a 2 r 2 



1 : 



r 



(p 2 -a 2 ) + X -(p 2 -a 2 ) 2 + a 2 z 2 , 



(34) 



(35) 



with p = \J x 2 + y 2 + z 2 . The difference here is the addi- 
tion of angular momentum. We test two cases, a = 0.5M 
and a = 0.9M. Fig. (7) presents a Schwarzschild (a = 
0M) case, together with the a = 0.5M and a = 0.9M 
cases. The solid line is the 8 = ir/2 slice and the dashed 
line is the (j> = tt slice. We find the expected result, that 
the deformation in the d> = tt slice increases with a. 






FIG. 7. The three plots correspond to the location of the 
apparent horizons for black holes with three different values 
of angular momentum. The units of the graph are M, the 
solid line is the 8 = -k/2 slice and the dashed line is the 4> = -k 
slice. 
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The radius of the event horizon is given by 

r+=M + \/M 2 - a 2 . (36) 



a critical separation 1.53(5)M. The apparent horizon at 
the critical separation of 1.535M is shown in fig. (8) using 
the level flow code with 33 2 grid points. 



The solution to eqn.(36) for a = 0.5M is r+ = 1.87M and 
the numerical solution we obtain for the horizon radius 
is r num = 1.87M, with an error of 0.06%. In the a = 
0.9M case, r + = 1.44M and r num = 1.46M, with a 
1.39% error. The area of the horizon for each case can 
be calculated using 

A = 4 7 r(r 2 _ + a 2 ) (37) 

(generalizing eqn.(30)), and compared to a numerical 
eqn.(31). computation using eqn.(31). In the a = 0.5M 
case, eqn.(30) gives A = 46.89M 2 , numerically we ob- 
tain A num = 46.88M 2 , resulting in a 0.21% error. In 
the a = 0.9M case, A = 36.09M 2 and A num = 36.39M 2 
with a 0.83% error. The errors will decrease as k is driven 
closer to zero. 



C. Brill-Lindquist Data 

In this section, we study a binary black hole system 
using Brill-Lindquist data [21]. These data are useful 
to us for two reasons: We can verify previous results of 
the critical separation, and study an example of how the 
tracker works in finding multiple apparent horizons. The 
3-metric is time symmetric, K a b = 0, and is conformally 
flat: 

9ab = <p 4 Vab (38) 

where 

N 

^=1 + E^ (39) 

i—l 1 

and TV is the number of holes (here TV = 2), Mj is the 
mass of the ith black hole and the are the radial dis- 
tances from the centers of the black holes 

We use isotropic coordinates to express the metric as 

ds 2 = 4> 4 (dr 2 + r 2 d8 2 + r 2 sin 2 8d<fi ) (40) 

with 



n = ^r 2 + d 2 - 2d ir cos 0, (41) 

where di is the distance between the holes and the center 
of the coordinate system. When they are far apart, each 
hole has an apparent horizon of radius M/2 in these co- 
ordinates. The area of each of the holes when they are 
well separated is 167rM 2 . 

The limiting separation of the holes between single and 
double horizons was found by Brill and Lindquist [21] to 
be 1.56M, Cadez 1.534M±0.002M [22], [1.5M, 1.6M] by 
Alcubierre et al. [23], and 1.535M by Huq [13]. We found 
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d= 1.535M 
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-1 1 

FIG. 8. Separation of 1.535M with a 33 2 grid. The area 
was determined numerically to be 184. 16M 2 

The horizon found for a separation of d = 1.5M which 
is less than the critical separation, is shown in fig. (9). 
Fig. (10) is a plot of the ^-norm of the maximum of k(9) 
for the separation d = 1.5M at each iteration plotted 
versus the number of A-steps. This is one of the checks 
in the level flow code to ensure that the apparent horizon 
equation is still being solved. We expect the expansion 
to continue to decrease if we have started outside the 
apparent horizon and are flowing inward. As we will see 
in fig. (14), the expansion increases as fission occurs in a 
data set with separated holes. 

, i , , , , | , , , , | , , 

, _ d= 1.5M 




-1 - 



-1 1 

FIG. 9. Separation of 1.5M with a 33 2 grid. The area is 
185. 41M 2 . 
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FIG. 10. The absolute value of the maximum of expansion, 
k, per iteration, A every 20th step. The kinks at A = 1000 
and A = 2000 are from restarting the code with a different 
A-step. 

As we increase the separation between the two holes to 
a separation greater than the critical separation, we can 
test the apparent horizon tracker in the case of multiple 
apparent horizons. We demonstrate with a separation 
of d = 2.0M. The initial surface flows to the point of 
fissioning where the topology of the surface changes from 
a one surface into two surfaces. Fig. (11) is a plot of 
the initial surface that begins the flow. The level set 
found during this flow is depicted in fig. (12). Each of the 
surfaces in fig. (12) has a constant expansion, k = c n and 
was used to indicate a topology change in the test surface. 
The values for the expansion are C\ = 0.14, c<i = 0.12, 
C3 = 0.1, and C4 = 0.08. The last single surface just 
before the topology change is not a surface in the level 
set; it is plotted in fig. (13). At this point the tracker 
begins to flow two surfaces. 




-2-1012 
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FIG. 11. The starting surface of the level flow method for 
a separation of d=2.0M 
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FIG. 12. The level set of surface for the d=2.0M case. Each 
surface has a constant k = c n at each point. In this case 
the values for c n are c\ = 0.14, C2 = 0.12, C3 = 0.1, and 
C4 = 0.08. The level set is used to indicate the change in 
topology associated with multiple surfaces. 

In contrast to a separation of d = 1.5M where there 
is no fission, here as fissioning becomes imminent, the k 
begins to increase. Fig. (14) is a plot of the absolute value 
of the maximum across the surface of the expansion, k, 
versus A up to the point of fission. The increase in the 
expansion is one of the signals of imminent fission. As the 
algorithm tries to find a surface with k = everywhere, 
it is driven into two surfaces. Once the new surfaces are 
found, the maximum of the expansion begins a monotonic 
decrease as in fig. (10). 



2.0M 
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FIG. 13. The single surface is about to fission into two 
surfaces. 
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FIG. 14. The absolute value of the maximum of expansion, 
k, per iteration is plotted. The increase in the expansion is 
caused by imminent fission. 

The exaggerated peanut shape in fig. (15) and fig. (16) 
is taken for the same A- value as fig. (13). 

Once the fissioning is detected by the code, it auto- 
matically begins flowing two new surfaces of the same 
resolution as the parent surface. The series of snapshots 
shown in fig. (15) and fig. (16) is a subset of the set of 
surfaces found by the apparent horizon tracker as it fol- 
lows the fission of the trial surface into two surfaces. The 
tracker starts with a spherical starting surface that de- 
forms along the gradient field. 



FIG. 15. This series of snapshots depicts the flow of an 
initial surface until its fission for the binary Brill-Lindquist 
black holes separated by 2M. The lower left plot is a first 
try at determining the final two surfaces. The cusps are due 
to a typical drawback associated with using points to define 
the flowing surface. The points crowd together in regions of 
greater flow. The next snapshot, on the lower right, shows the 
code's automatic correction; and shows the apparent horizons 
of the binary Brill-Lindquist data to an accuracy of 10~ 4 . 

As the points defining the surface flow, the distance 
between the points can become too small for the finite 
difference scheme at that resolution. Redistribution of 
the points on the surface is taken care of automatically 
by updating the center and radius. 




-2-1 1 

K[Ml 
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FIG. 16. The series of pictures shown in fig. (15) are placed 
in one plot. The outer surface is the initial guess, the "peanut 
surface" is the surface that is found indicating the need to 
search for two surfaces, and the inner approximate spheres 
result from locating the apparent horizons for Brill-Lindquist 
data. 
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VII. APPARENT HORIZONS IN A GRAZING 
COLLISION 

As stated in the Introduction, one of the main moti- 
vators of this work is to have an apparent horizon finder 
that can locate disjoint horizons during the evolution. 
This entails 1) finding the horizons without a good ini- 
tial guess, and 2) detecting the topology change from two 
disjoint horizons to one horizon. To demonstrate the level 
flow's ability to carry out 1) and 2), we report the results 
of apparent horizon tracking in the particular case of a 
short evolution of two spinning, Kerr-Schild black holes 
using the Binary Black Hole Grand Challenge Alliance 
Cauchy code [24]. A future paper [25] will discuss the 
details of the evolution. 
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The evolution is free, i.e. the momentum and Hamil- 
tonian constraint equations are only used as checks dur- 
ing the evolution. Since we cannot hold infinity on the 
grid in this formalism, we must specify outer boundary 
conditions for the dynamic variables, gij and Kij. For 
the following work, we specify analytic outer boundary 
conditions with blending [26] between the analytic and 
numeric regions. 

To specify initial data for two spinning, boosted black 
holes we use superposed Kerr-Schild black holes. We 
chose a Kerr-Schild metric [20] for two reasons: 1) The 
metric is well defined at the event horizon, and 2) The 
metric is Lorentz form-invariant in a simple sense, under 
boosts (v > 0). The superposed data are constructed 
in an approximate manner by a conformal method based 
on the superposition of two isolated, boosted Kerr black 
holes. 

The initial data follows from Matzner et al. [27] and 
was first implemented by Correll [28]. The data is the 
superposition of two, isolated, boosted Kerr-Black holes 
with individual metrics given by eqn.(23). The resulting 
superposed metric is: 

9ij = (i)9ij + (2)9ij ~ Vij (42) 

with the " symbol indicating a quantity conformally 
related to the physical metric, <?y = (f> 4 gij. 

(\)9ij = Vij + (i) H ( r i)(i)k(i)lj and (43) 
(2)9ij = Vij + (2)H(r 2 )( 2 )li(2)lj (44) 

are the the isolated Kerr-Schild metric forms with l t and 
H corresponding to the single black holes. The two holes 
have comparable masses, Mi ~ M 2 , coordinate separa- 
tion ri2, and velocities v\ and v 2 assigned to them. For 
the argument of H and lj , we use 

r\ 2 = (x — x\) l (x — x\) 3 8ij and (45) 
r 2 2 = (x-x 2 ) i (x-x 2 yS ij (46) 

with x\ and x-j the coordinate positions of the holes on 
the initial slice. 

The extrinsic curvatures of the two isolated black holes 
are added to obtain a trial K a b for the binary black hole 
system given as: 

(o)Kij = (\)Kij + ( 2 )Kij. (47) 

The subscript indicates that this is an approximation 
to the true extrinsic curvature of the binary black hole 
spacetime. and ( 2 )Kij are the individual extrin- 

sic curvatures associated with the isolated Kerr-Schild 
metric and their indices are raised and lowered by their 
individual metrics, eqn.(43). 

For the data we describe here (holes center initially 
separated by a coordinate distance exceeding 10M where 
M is the mass of one of the holes), we expect that an 
initial value solution will be £s 10% in error on the domain 



outside of the excision volume. See further discussion in 
[29]. We set the lapse function to: 

a = ai + a 2 - 1, 

and the shift vector to: 

F = p[+&. 

The run presented in this paper has a grid 81 x 81 x 
81 in Cartesian coordinates (x,y,z) with a domain of 
(±10M, ±10M, ±10M) resulting in a spatial resolution 
of M/4. The data represent two black holes in a grazing 
collision. The holes are set initially at (5M, 1M, 0M) 
and (— 5M, — 1M, 0M) in Cartesian coordinates with a 
boost speed of ±0.5x toward one another and each has 
an angular momentum per unit mass of a = 0.5M in 
the (— )z-direction. Fig. (17) is the initial configuration of 
this run; note that a naive sum of the spin and the orbit 
angular momentum yields zero for this configuration. 

We post-process the data obtained from the evolu- 
tion. For the purposes of this paper, we track the appar- 
ent horizons at three specific times during the evolution, 
namely t = 0M, 2.8M, and 3.4M. At 0M, the apparent 
horizons of the initial data are found, at 2.8M two dis- 
joint apparent horizons are found; and finally, at 3.4M a 
single apparent horizon is found. For the horizons shown 
here, the level flow method used a sphere of radius 8M to 
initialize each run. Fig. 18 is a plot of the horizons with 
time going up the page. The lowest plot is of t = 0M, 
with each horizon being a sphere centered at coordinates 
(±5M, ±1M, 0M). The middle plot shows the horizons 
at a later time, t = 2.8M. Here the deviation in shape as 
the horizons accelerate towards each other is seen. The 
final plot at the top of fig. 18 is the first single apparent 
horizon that envelops both black holes at t = 3.4M. 

The areas for the apparent horizons at t = 0M are 
A = 43. 6M 2 for each hole. At t = 2.8M, the horizons 
have deviated from a spherical shape and the areas for 
each hole are A = 44. 2M 2 , giving a measure of the ac- 
curacy to which we can maintain their areas constant. 
The area of the merged apparent horizon at t = 3AM is 
Amerged = 184M 2 . According to the black hole area the- 
orem of Hawking and Ellis [30], the area of the merged 
event horizon must equal at least the combined area of 
the individual event horizons. Although no strong state- 
ments can be made about the area of an apparent hori- 
zon, we do find A merge( i > A\ + A 2 in a consistent man- 
ner. We can further surmise that the final maximum area 
we could expect based on the initial configuration should 
be that of a Schwarzschild black hole of mass 2M, giving 
an area of approximately 201M 2 . In some sense the area 
predicted by the Schwarzschild case is an upper bound. 
We see a 8.5% deviation from that "idealized" case. In 
view of this upper limit, 8.5% may be an indication of 
the greatest amount of gravitational radiation up to the 
time of merger (t = 3.4M) given our approximate initial 
data, gauge condition, and boundaries. 
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VIII. CONCLUSION 



a=-0.5z 



y[M] 



v=0.5x 



x[M] 



FIG. 17. The configuration of the initial data for the graz- 
ing collision. The initial mask position is indicated by the 
"circle" centered on each hole. The angular momentum per 
unit mass and velocity of each of the holes is also represented. 
The total angular momentum (spin plus orbit) of the initial 
configuration is zero. 




Apparent horizon location and tracking constitute an 
important part of numerical evolutions of black hole 
spacetimes using excision techniques. We have demon- 
strated a method for finding apparent horizons in situa- 
tions where the location of the apparent horizon may not 
be known; hence a good initial guess for the finder may 
not be possible. The method we have discussed works 
with generic 3-metric and extrinsic curvature data and 
with an arbitrary initial starting surface. Furthermore, 
the method is capable of detecting a topology change 
as the finder flows towards the apparent horizon. This 
ability is important for situations where there are mul- 
tiple apparent horizons in the data. This allows us to 
locate apparent horizons in binary black hole evolutions 
without knowing where the apparent horizons are; and it 
allows us to locate the first single apparent horizon that 
forms at the merger of two black holes. The level flow 
method is successful at locating the apparent horizons in 
generic spacetimes as demonstrated by the Schwarzschild 
and Kerr data. It also found multiple apparent hori- 
zons starting from a single starting surface as demon- 
strated with the Brill-Lindquist data. Most importantly, 
the level flow method has been successful at identifying 
apparent horizons in a binary black hole evolution involv- 
ing two Kerr-Schild black holes. Beginning with a single 
guess surface, two discrete apparent horizons were found 
at early times, and the later single merged horizon was 
found. 

One of the drawbacks of this method currently is its 
slow convergence property due to the parabolic nature 
of the equation solved. This, however does not pose a 
problem since the level-flow method can be used in con- 
junction with other methods which may be more efficient 
given a good initial guess. The level-flow method has the 
definite advantage of being capable of finding multiple 
surfaces in the data. It can be used to get extremely 
good initial guesses for other methods that converge more 
quickly close to the solution. We are currently using the 
level flow method in this manner in numerical evolutions 
of black hole collisions. 




FIG. 18. The apparent horizons are plotted with the evo- 
lution time increasing up the page. The times plotted are 
0.0M, 2.8M, and 3AM. 
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